Non-equilibrium time-dependent solution to discrete choice with social interactions

We solve the binary decision model of Brock and Durlauf (2001) in time using a method reliant on the resolvent of the master operator of the stochastic process. Our solution is valid when not at equilibrium and can be used to exemplify path-dependent behaviours of the binary decision model. The solution is computationally fast and is indistinguishable from Monte Carlo simulation. Well-known metastable effects are observed in regions of the model’s parameter space where agent rationality is above a critical value, and we calculate the time scale at which equilibrium is reached using a highly accurate method based on first passage time theory. In addition to considering selfish agents, who only care to maximise their own utility, we consider altruistic agents who make decisions on the basis of maximising global utility. Curiously, we find that although altruistic agents coalesce more strongly on a particular decision, thereby increasing their utility in the short-term, they are also more prone to being subject to non-optimal metastable regimes as compared to selfish agents. The method used for this solution can be easily extended to other binary decision models, including Kirman’s model of ant recruitment Kirman (1993), and under reinterpretation also provides a time-dependent solution to the mean-field Ising model. Finally, we use our time-dependent solution to construct a likelihood function that can be used on non-equilibrium data for model calibration. This is a rare finding, since often calibration in economic agent based models must be done without an explicit likelihood function. From simulated data, we show that even with a well-defined likelihood function, model calibration is difficult unless one has access to data representative of the underlying model.


Introduction
It has been 20 years since the publication of discrete choice with social interactions [1], which at the time of writing has over two thousand citations. The success of the publication has spawned a multitude of related publications, all with an interest in modelling how variably rational economic agents make decisions under exogenous influence in a system with collective endogenous interactions. More explicitly, the model of Brock and Durlauf considers a system of agents where each agent makes a decision left or right, where their decision is influenced by global influences (affecting all agents) and collective effects relating to a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 However, for all the applause of the model of Brock and Durlauf, the solution that they provide for the binary choice model applies only at equilibrium, whereas it is now widely known that the time to reach equilibrium, for many socio-economic systems of interest, does not occur on an economically relevant time scale [3,25]. Hence, to understand binary decision phenomena more deeply one must consider the whole time-evolution of the dynamics from an initial condition to the equilibrium. In particular, model calibration is likely to be greatly perturbed using a likelihood function based on the equilibrium solution for data trajectories that are not a sample of the equilibrium distribution. In this paper we solve for the probability distribution of the mean-field binary decision model of Brock and Durlauf in time, and use our solution to investigate metastability, and the requirements on economic data necessary for model calibration. The method used for our solution also solves Kirman's ant model for the case of discrete numbers of agents, complementary to the study of [11], where a time-dependent solution is provided for a continuous number of agents (i.e., the large agent population limit). This analytic solution allows for faster model calibration compared to simulation based approaches, as well as the analytical construction of a likelihood function.
The paper is structured as follows. In Section 2 we introduce the mean-field binary decision model of Brock and Durlauf starting from the RFIM of Bouchaud [3]. In this way the reader can see how the mean-field description comes from applying various assumptions to the RFIM. Section 2.1 describes the situation of a system of selfish agents who make decisions such that they maximise their utility (if they are rational). In Section 2.2 we extend this analysis to a system of altruistic agents, where now agents are concerned with optimising the global utility. Then in Section 2.3 we solve the mean-field binary decision model in time. Throughout Section 3 we explore the solution that we have found the metastable phenomena it possesses (known as a lock-in in the economics community). In the process, using a first passage time analysis, we outline a method to very accurately determine the time scales of these metastable states. In Section 4 we use our time-dependent solution to construct a likelihood function that can be used for model calibration, and assess the conditions necessary on simulated data to provide accurate calibration. We then compare the results from our model to the standard tenets of neoclassical economics in Section 5, before concluding our paper in Section 6.

Selfish agents
Consider a system of N economic agents where each agent i can choose between two distinct decisions, left denoted by S i (t) = −1, and right denoted by S i (t) = 1. Left or right could correspond to an infinite number of binary choice questions, for example, to vote either Democrat or Republican, to buy one stock or another, or whether one is participating in the latest fashion trend (in this case yes or no). Depending on the question one is considering, the model parameters making up a binary decision model are likely to be very different. For example, a binary decision model of an American election is unlikely to have the same model parameters as one describing consumers choosing between two brands of similar cereal: in the case of a cereal, the exogenous influence of advertising is more important than endogenous effects, whereas the collective effects are likely more important in an election model. The attribute S i (t) denotes the decision currently made by agent i at a time t, with the state of the entire system at time t being given by the set of choices made by all agents at that time, denoted S ¼ fS 1 ; S 2 ; . . . ; S N g. The RFIM introduced by Bouchaud in [3] takes into account three main factors, all of which contribute to the influence on agent i, I i (t), which are (using naming conventions inspired by [3]): 1. The personal inclination (or predisposition) of the agent to favour one decision over the other. The contribution to I i (t) from personal inclination is f i 2 [−1, 1], where f i < 0 means a predisposition to left and f i > 0 means a predisposition to right. In principle, these heterogeneities could be time-dependent, but since they are difficult to empirically quantify one generally assumes they are constant.
2. The zeitgeist, a global influence equally affecting all agents, such as news reporting or stock trends. Its contribution to I i (t) is given by F(t) 2 [−1, 1], and we assume its time-dependence directly since it is often global effects (e.g., stock market decline or warming global temperature) that change in time.
3. Agent-to-agent interactions, attractive/repulsive influences for the agents to attract/repel others from their same view. The strengths of these interactions are stored in the matrix J, where J ij � [J] ij is the contribution of agent j to I i (t) such that the total contribution of all agents in neighbourhood ν i to the influence of agent i is P j2n i J ij S j ðtÞ.
In summary, each agent is subject to an influence, I i (t), which is the sum of these contributions, The importance of the influence is that, where agents are rational, they will aspire to make a decision (left or right) that agrees with the influence they are subjected to. We explore the mechanisms through which agents update their decisions below, based on the introduction of an agent utility function.
In the mean-field case several simplifications can be made. The main assumption is that all agents are equally connected to each other with strength J meaning that each agent only responds to the average global opinion of all agents m = ∑ i S i /N. As has been widely documented, the mean-field approximation is qualitatively a good approximation for nearest neighbour interactions where N � 1 and the number of dimensions d � 4 or in situations where agents are connected to many other agents [26]. We additionally rescale J by the number of agents J ! J/N, such that the agent-to-agent interaction contribution to I i ðS; tÞ is intensive. Finally, we assume that each agent has the same personal inclination f, a factor that is absorbed into the definition of the F(t) in Eq (1). Note that even in the absence of personal preference the agents are still heterogeneous in the sense that they make their own decisions at different times. Within this set of approximations the information is a shared function amongst all agents and is given by, where we reinterpret the definition of m(n(t)) = (2n(t) − N)/N, where n(t) is the number of agents whose decision at t is right (with N − n(t) deciding left). For brevity, we will often drop the time-dependence of n on t, although it is always implicit. The function m(n), as well as being the average global opinion, is the order parameter of the system and takes discrete values in [−1, 1] for finite N: m(N) = 1 if all agents take the right decision, m(N/2) = 0 if equal numbers decide right and left, and m(0) = −1 if all agents decide left. Hence, m(n) appropriately characterises the state of the system for any value of N.
In the mean-field context one can interpret the meanings of F(t), J and m(n) in many ways. One way, introduced in [3], is to interpret F(t) as the cost of installing a new heating system and Jm(n(t)) as the expected cost decrease as technical improvements cause more people to switch to the dominant technology. A similar interpretation will be used later on in Section 3.1 to discuss the phenomena of technology lock-ins. In many settings, mean-field approximations are reasonable for binary choice decisions: many decisions one makes on a day-to-day basis are influenced by many peers acting around an individual, for example, deciding whether to invest in cryptocurrency [27], or the language one uses in online reviews [28]. As noted by Kirman [6], sociologists have long observed empirically that relational networks are likely to be much more connected than one might imagine. However, in many cases network effects are very relevant, for example, children brought up in religious households are much more likely to be religious themselves [29], in which case the mean-field assumption is a poor one. We illustrate the concept of the mean-field assumption in Fig 1(a).
We now introduce the utility of agent i as U i ðn; tÞ � S i ðtÞI i ðS; tÞ ¼ S i ðtÞIðn; tÞ, a quantity that a rational agent would want to maximise over any realisation of the model. We denote the situation wherein agents only care to maximise their own utility as a system of selfish agents and generalise for more altruistic agents that are motivated by changes in the global utility in the following section. Now consider that the rate at which an agent changes their mind is Fig 1. (a) Illustration showing the effect of mean-field theory on a network of interconnected agents that influence each others decisions. In the meanfield case each agent is influenced by the same information I(n, t) which they all contribute towards. (b) Diagram showing the 1D process that is equivalent to the mean-field case of the binary decision model. The number in each circle denotes the number of right deciding agents, and the expressions above and below the arrows denote the associated probabilities of transitioning between each configuration. (c) Plots showing the correspondence between our analytic solution in Eq (16) and Monte Carlo simulations provided by the SSA. The parameters for this time-evolution are F = 0, J = 10, α = 0, β = 1, γ = 1, N = 100 and P(n, 0) = δ n,N/2 . The probability distribution from the SSA is calculated from 2.5 × 10 3 trajectories. As t ! 1 we see the emergence of a single steady-state consisting of two equal modes of height �0.5 at m(0) = −1 and m(N) = 1. In the limit N ! 1 this steady-state bimodality corresponds to symmetry breaking behaviour [30]. Note that the SSA simulations in the solid lines show random perturbations due to the stochastic nature of the simulations, whereas the bars showing Eq (16)

PLOS ONE
dependent on the change in their own utility if they do change their mind. For agent i, whose utility at time t is U(n, t), changing their decision from S i ! −S i at time t would result in the utility change, Note that 2J/N is a self-interaction term which comes from considering the change in m(n) as S i ! −S i , a contribution that equally affects the rates of transition from −1 ! 1 and 1 ! −1, and that is negligible for N � 1. In the case where the change in the agent's decision is not assumed to affect ΔU i (n, t) the term 2J/N can be ignored. If an agent had made decision S i (t) = −1 (or S i (t) = 1) and the influence on them was I(n, t) > 0 (or I(n, t) < 0) then in changing their decision at t they increase their utility; alternatively, if an agent had made decision S i (t) = −1 (or S i (t) = 1) and the influence on them was I(n, t) < 0 (or I(n, t) > 0) then in changing their decision at t they decrease their utility. Hence, a rational agent will (on average) make decisions such that their choice S i (t) agrees with the sign I(n, t), at a rate (determined below) that is a monotonically increasing function of ΔU i (n, t). More simply, the larger the change in the agent's utility ΔU i (n, t) upon a decision change, the larger the rate at which an agent changes their decision.
In this paper we choose decision rules based on detailed balance considerations [30]. This means that we choose the rates with which agents update their decisions based on the assumption that as t ! 1 the Boltzmann equilibrium distribution will be attained, i.e., P eq (S i ) / exp (βU i (n)) where β is the rationality of agent i (discussed below). Enforcing detailed balance means the transition rates must satisfy, where W n (S i (t) ! −S i (t)) is the transition rate for agent i to change their decision from S i ! −S i at time t given there are already n agents deciding right [30]. That W n (S i (t) ! −S i (t)) is a transition rate means that in a time interval [t, t + Δt) an agent will change their decision with probability W n (S i (t) ! −S i (t))Δt. There are several conventional ways to prescribe the transition rates from Eq (4). A common way is to choose the Arrhenius form [31] where the transition probabilities become W n (S i (t) ! −S i (t)) = γ exp(βU i (n − S i , t)) and W nÀ S i ðÀ S i ðtÞ ! S i ðtÞÞ ¼ g expðbU i ðn; tÞÞ, where γ is a time scale parameter. An alternate form, and the form we use in this paper, are the Glauber/logit transition rates, which are given by [3,32,33]: The choice of logit transition rates is motivated by [33], which derives Eq (5) using the maximum entropy principle. This assumes that agents make their choices based on a compromise between short-term utility gain and a desire to sample other available choices. We note that the Arrhenius transition rates do not subscribe to this interpretation, since as is clear from their form, the transition rates are only dependent on the utility of the agent upon the change in decision; i.e., the agents do not consider the change in decision based on a knowledge of all decisions they could take, only the one they will take. In reality the proper form of the transition rates likely varies depending on the problem at hand. However, it seems likely that in the case of a binary decision an agent would need to first explore a range of decisions before they know the one that maximises their utility, hence we use Eq (5) for the rest of the paper. In other models the logit-like behaviour arises from other considerations, for example their presence in the future technology transformation (FTT) models described in [34] is attributed to probabilistic nature of the cost of two competing technologies (see Fig 1 of [35]). One also observes that the dependence of the propensity function W n on the change in utility is identical to that seen in classical choice theory [36]; this is by design, and the propensity of an agent to change their decision is proportional to the probability that the agent would indeed choose that option in choice theory. Eq (5) allows us to explore the agent rationality β, which is a direct analogue of the inverse temperature commonly seen in equilibrium thermodynamics and statistical physics [30]. As β ! 0 the agents become completely irrational and changes of decision occur at the same rate γ/2 regardless of an agent's utility change. Conversely, as β ! 1 each agent is completely rational and Eq (5) becomes, ( Hence, completely rational agents always make decisions that agree with the influence upon them. For intermediate values of β we obtain the typical 'S-shaped' adoption curves with respect to ΔU i (t), a common feature of many macroeconomic models including heterogeneous agents [34,37,38]. We note that since β is a constant we are implicitly assuming that each agent has the same level of rationality. Although in reality it may not be the case, for most situations we assume it to be a good approximation.
In the model described above 'completely rational' agents correspond to agents that on their next change of decision, only choose the option that maximises their utility. So in a sense, even completely rational agents have limited foresight since they make decisions only based on the current state of the system, and not based on which decision in the long-run will optimise their utility-which would be the decision agreeing with the sign of F. Clearly, even these hyper-rational agents do not agree with the 'perfect rationality' of agents observed in neoclassical economics (where agents even optimise for future conditions [39]). The agents in the binary decision model more closely correspond to Keynesian agents with 'fundamental uncertainty' in economies driven by short-term optimism (in our case short-term utility gain) [40].
In contrast to the above, Kirman's model of ant rationality [10] assumes a different form of propensities than the logit form seen in Eq (5). Interpreting S i (t) = ±1 now as two different sources of food and n as the number of ants at the right-hand (+1) food source, ants switch food source due to two influences: the first is due to random switching with probability � per unit time, the second due to recruitment by another ant with probability μ per unit time. The former is akin to an effective reaction with first-order kinetics, whereas the latter is akin to a reaction with second-order kinetics-i.e., the chance meeting of two ants is modelled as the probability of two molecules of the same species colliding. Following the laws of mass-action kinetics [14,31,41], the propensities for Kirman's ant model are then proportional to W n ðÀ 1 ! þ1Þ ¼ ðN À nÞ� þ mðNÀ nÞn NÀ 1 and W n ðþ1 ! À 1Þ ¼ n� þ mðNÀ nÞn NÀ 1 , as seen in [11]. We will see below in Section 2.3 that our method for solving the model of Brock and Durlauf in time also easily extends to Kirman's model, although we focus on the former model for the rest of the paper.

Altruistic agents
We now generalise our mean-field model of economic agents such that they can take into account changes in the global utility of all agents. We refer to this as agent altruism since in this case agents do not only think about maximising their own utility but that of the entire population [42]. In the context of binary decisions, having agents with altruistic attributes is somewhat realistic. Take the example we made above relating F(t) to the cost of a heating system and Jm(n) as the reduction in cost given more agents are using the technology. In this example altruistic agents could further benefit the utility of all individuals by more quickly arriving at a state where all agents are in agreement with one another. We define the global utility simply as the sum of the utilities of all agents, Due to the mean-field interactions between the agents, a change in decision in one of the agents will change the utility of the individual agent by a different amount than the global utility. For example, say we have S 1 ¼ fS 1 ; . . . ; S j ; . . . ; S N g and at time t the system evolves to , unless there are no interactions between the agents. We consider A(n − S i , t) − A(n, t) explicitly below. Now consider a system of entirely altruistic agents, whose decisions do not depend on their individual utilities but only on A(t). As we did before for the selfish agents, we want to calculate the change in the A(t) if agent i changes their decision from (6), one finds that this is given by, Hence, the change of the global utility is equal to the change in the individual utility of the agent minus a non-negligible term that accounts for the change in the utility of the rest of the population given agent i has changed their decision. Note that the m(n) − S i /N term comes from considering the sum over all agents aside from agent i, The extra term in Eq (7) has an intuitive meaning, whose presence is explained by two mechanisms: (i) if S i (the original choice made by the agent) has the same sign as m(n) then ΔA i < ΔU i since the agent has decided to go against the average opinion of all other agents encoded in m(n); (ii) if S i has the opposite sign to m(n) then ΔA i > ΔU i since the agent has decided to go with the average opinion, having previously been against it. Following previous work [42] we now introduce the gain, which gives a generalised form of utility change for agents that can be either entirely selfish, altruistic, or else somewhere in between. The gain for agent i in changing their decision from S i ! −S i is defined by, where α is a parameter such that α = 0 corresponds to a system of selfish agents, α = 1 a system of altruistic agents, and 0 < α < 1 somewhere in between. We can finally write our most general transition rates including the effects of agent altruism, which is essentially identical to the decision rule of Eq. (1) in [42]. Curiously, and possibly somewhat intuitively, including the effects of altruism in our mean-field description acts to scale the agent-to-agent interaction strength by 1 + α. As we explore in Section 3, this scaling of the endogenous strength can have qualitative as well as quantitative effects on the model. It seems the mean-field approximation is responsible for the simplicity of this scaling, since in the general case of the model introduced in Section 2 the difference in utility change between altruistic agents and selfish ones is found to be the less simplistic and given by Note that for the case of generalised mean-field altruistic agents one can define a Hamiltonian HðSÞ, a function of the state of the system S such that, One finds this function is given by, where we remind the reader that m(n) = (2n − N)/N = ∑ i S i /N. The existence of HðSÞ ensures that decision rules based on Eq (9) satisfy detailed balance and that as t ! 1 the equilibrium probability of each state is P eq ðSÞ / expðbHðSÞÞ for all 0 � α � 1.

Master equation and analytical time-dependent solution
We can now proceed to solve the mean-field model analytically. First we map the model to a stochastic birth-death reaction scheme given by, where r(n) and l(n) are the transition rates at which each agent changes their decision from left to right and right to left respectively given there are n right voting agents, and L and R are symbols representing agents who have decided left and right respectively. In this context one can state that a decision change of an agent is akin to a gain/loss of particles, where each particle denotes a right deciding agent [31]. Importantly, from this birth-death mapping, where F(t) ! F is a time-independent constant, it is clear that in the mean-field limit the number of right deciding agents n completely determines the state of the system. Hence, for the moment we set F (t) ! F, and discuss the case of time-dependent F(t) later on. The rates r(n) and l(n) are defined by the decision rules we derived above in Eq (9) and are explicitly given by, Note that the choice of the overall propensities (N − n)r(n) and nl(n) comes from mass-action kinetics: in short, the rate of gaining a new right deciding agent is proportional to the number N − n of left deciding agents in the system [14,31].
One can easily check that by writing the rate equation for the evolution of Eq (12), in the large N limit, as t ! 1 one recovers the classic mean-field result of m = tanh(β(F + J(α + 1) m)). However, we can do better than the rate equation describing the steady-state value of m, and can in fact solve for the probability of having n right deciding agents at a time t, Pðn; tÞ. We first write the master equation corresponding to the process in (12), @ t Pðn; tÞ ¼ ½ðN À ðn À 1ÞÞrðn À 1Þ�Pðn À 1; tÞ þ½ðn þ 1Þlðn þ 1Þ�Pðn þ 1; tÞ À ½ðN À nÞrðnÞ þ nlðnÞ�Pðn; tÞ: The master equation, otherwise known as the Kolmogorov forward equation, is a set of firstorder differential equations describing the time evolution of being in discrete state n at time t.
Eq (14) is a 1D master equation since it describing the evolution of only a single stochastic variable n. We outline the derivation of a 1D master equation from first principles in Sec. 1 in S1 File. At steady-state, where @ t Pðn; t ! 1Þ ¼ 0, 1D master equations are very well understood, and their properties are discussed in many textbooks [30,31,43]. Solving them in time is a much trickier problem. However, it is possible to solve Eq (14) and we do so using the non-standard and useful method of [44]. For readers that are less mathematically inclined we suggest that not too much time is spent on the derivation below; it is simply most important to realise that the mean-field binary decision model can be solved in time, and that its analytical solution is much quicker to compute than simulation based methods. For consistency with [44] publication we introduce the following notation: a n = (N − (n − 1))r(n − 1) and b n = (n + 1)l(n + 1), and we show the transitions between the microstates of the model using these propensities in Fig 1(b). Then, Eq (14) can be re-written as @ tP ðtÞ ¼ A �P, wherePðtÞ is the column vector ðPð0; tÞ; Pð1; tÞ; . . . ; PðN; tÞÞ and A is the master operator, a (N + 1) × (N + 1) tridiagonal matrix, given by, We denote the eigenvalues of the matrix A as λ i for i = 1, 2, . . ., N + 1, and from Perron-Frobenius theorem can state that the largest eigenvalue is [45]. Note that the eigenvector corresponding to λ 1 = 0 is the steady-state probability distribution. Following some complex analysis on the resolvent of the master operator and judicious use of Cauchy's residue theorem in [44] one then arrives at the solution of Eq (14) from an initial condition Pðn; 0Þ ¼ d n;n 0 , explicitly, where the orthogonal polynomials p n and q n are recursively defined via, Eqs (16)- (20) constitute the analytical time-dependent solution for the mean-field Ising-Weidlich binary decision model with variably selfish-altruistic agents [46]. Note that unlike Pðn; tjn 0 ; 0Þ, P(m(n), t|m(n 0 ), 0) is not a probability density but instead the raw probability of the system having a certain value of m(n) at time t, i.e., ∑ n P(m(n)) = 1 and ∑ n P(m(n))Δm 6 ¼ 1 (where Δm = 1/N). For details regarding the derivation of this result we refer the reader to [44]. In the most general case of usage of the model one needs to determine the eigenvalues of A computationally, which we do using the eigvals function in the Julia package Linear-Algebra [47]. However, since the eigenvectors are already implicit in the form of Eq (16) we do not need to evaluate these computationally, and hence the analytical solution we utilise can be computed much more quickly than a finite state projection approach requiring matrix exponentiation [48]. Note that the eigvals function that we use to compute the eigenvalues does not always return entirely real eigenvalues in every case (as one would physically expect) and often eigenvalues come in the form of a complex conjugate pairs. This is a common computational error, since we know theoretically that the eigenvalues should be real, and the resultant set of eigenvalues we obtain from eigvals is known as a pseudospectrum [49], which arises since the eigenvalues of these matrices are very sensitive to small perturbations. Importantly however, the usage of the pseudospectrum in Eqs (16)-(20) returns a normalised probability distribution that is indistinguishable from Monte Carlo simulation of the model. The Monte Carlo simulation method we employ is known as the stochastic simulation algorithm (SSA), and we detail it in Sec. 3 in S1 File. Briefly, the SSA is a continuous time method often used in stochastic chemical kinetics to simulate chemical reactions. In our case, it outputs individual realisations of the master Eq (14), which can then be replicated as an ensemble in order to construct the probability distribution, mean and variance of the process in time. In Fig 1(c) we show that our analytic time-dependent solution corresponds precisely to the distribution produced via the SSA for times near the initial condition through to near steady-state conditions. In Section 3.2 we approximately calculate λ 2 , the eigenvalue that determines the relaxation rate to the steady-state, using exact expressions for the switching times between lock-in states. Note that l À 1 2 gives the time scale to reach the equilibrium. We further note that an analytical procedure to calculate the eigenvalues of our problem in a perturbative sense may be possible following methods used in [50].
Eqs (16)- (20) assume that the initial condition is a fixed value of n 0 right deciding agents at t = 0. However, often the initial state is not precisely known but is itself a distribution QðnÞ � Pðn; 0Þ. Using the laws of probability one determines the time evolution of Pðn; tÞ with initial condition QðnÞ as, PðmðnÞ; tjQðmðnÞÞ; 0Þ ¼ Pðn; tjQðnÞ; 0Þ ¼ where QðmðnÞÞ ¼ QðnÞ. In the common case where each agent is initially assigned a decision at random with probability p 0 the number of initially right deciding agents is drawn from a binomial distribution, i.e., QðnÞ ¼ Binðn; N; p 0 Þ. The steady-state distribution reached as t ! 1, i.e., P s ðnÞ ¼ Pðn; t ! 1Þ, is widely known for birth-death processes with general rates and can be expressed as the following from Kirchhoff's theorem [51], where we define the empty product Q i j>i as being equal to 1. P s (n) is independent of the initial condition, even for systems with long-time metastability/lock-in effects. Note that for the case of time-dependent F(t) the solution above must be somewhat modified since a i and b i are no longer dependent only on n but also explicitly on t. This case of having a zeitgeist that changes in time has important policy implications, and although it is not relevant for our work here we outline its solution in Sec. 3 in S1 File.
Finally, we note that if one interprets the economic agents as magnetic spins, the zeitgeist instead as a magnetic field, and J as the mean-field exchange coupling (multiplied by the number of nearest neighbours), then our time-dependent solution above also provides the solution to the mean-field Ising model for a finite number of N magnetic spins. This model of a magnet, famed for its simple equilibrium solution among other things, is still the subject of recent work [52][53][54].

Exploration of the model
We can now use our analytic time-dependent solution to explore the model in different regimes of parameter space. Before delving into metastable behaviours, which will occupy us for the rest of this section, we detail two well known features of the model [1,3,30]. The first is that, for F = 0, the model exhibits a phase transition at a critical value of agent rationality β c = 1/J(1 + α). The origin of this critical value is explored in Sec. 4 in S1 File, but essentially for β > β c one finds a symmetry breaking behaviour wherein agents either mostly become left or right deciding at steady-state and P s (n) is bimodal; whereas for β � β c the steady-state distribution P s (n) is monomodal with a peak at m = 0. At β = β c we see the characteristic shape of a flattopped steady-state distribution centred at m = 0 such that any small deviation above β c causes bimodality. Secondly, for F 6 ¼ 0 one finds that although there may still exist two equilibrium modes of behaviour deterministically, in reality the decision of same sign to F is exponentially more favourable as t ! 1 (shown explicitly in Section 3.2). Economically, what this means is that agents coalesce on the same decision when they are sufficiently rational and when collective effects are deemed at least as important as exogenous effects. A recent review of pathdependency effects similar to lock-ins is given in [55].

Lock-ins
Let's explore a possible narrative of the model and suppose the two decisions left and right correspond to two competing technologies. F accounts for any exogenous effects on the agents, for example, the cost difference between the two technologies or the effects of changing social norms. Denoting C L and C R as the weight of these real and social costs of the left and right technologies, one can the break down F into its constituent parts, i.e., F = C L − C R . On the other hand, the interaction term J takes into account any endogenous changes arising from within the dynamics of the model, in particular: (1) learning effects as technology gets improved due to increased uptake [56], (2) any endogenous changes in price due to the increased uptake [56,57], and (3) establishment effects, accounting for the ease of investment regardless of the price since the infrastructure needed for the technology is already there [56]; i.e., J = J 1 + J 2 + J 3 .
Importantly, the foresight of the agents is limited only to how their decision change at the current time changes their own (or global) utility, and they do not know a priori which technology is the optimal one (i.e., they cannot distinguish between exogenous and endogenous influences). The state of maximal utility for the agents will clearly be that where all agents decide on the optimal technology (leading to the highest utility amongst all agents), which has the best combined trade-off between price and social appropriateness. However, it has been observed that this does not always happen, and path-dependency can lead to people making collective decisions that differ from the optimal one. This leads to so called lock-in effects [58], a type of positive reinforcement metastable phenomenon. As highlighted by [58], several types of lock-in are observed in the real-world, including technological, institutional and carbon lock-ins. Clearly, understanding lock-in effects is important for policy considerations. For example, where there are competing technologies with some being more carbon efficient than others, consumers will be more likely to invest in the cheapest most common technology with the most infrastructure than the one that is better for the environment [34].
In some regions of parameter space lock-ins occur in the model of Brock and Durlauf, where there are 2 possible lock-in states m � : the optimal one where sign(m � ) = sign(F), and the non-optimal one where sign(m � ) = −sign(F), which occur since the agent-to-agent interaction terms become dominant over the zeitgeist, i.e., switching to the alternate technology punishes the utility of an agent due to the collective effects. However, whether a lock-in is observed or not depends on several factors. First, one must be in the regime wherein agents are rational enough to begin coalescing on either of the two choices, i.e., β > β c . In this sense, β c is not only a property of the steady-state but also a requirement for certain types of transient phenomena. Second, F should not be of too great a magnitude such that only 1 equilibrium solution of m is present (see Sec. 4 in S1 File). Thirdly, the initial condition should not be entirely weighted towards the steady-state distribution one finds as t ! 1, i.e., since lock-ins are a path-dependent phenomenon there must be a viable opportunity for them to occur.
In many complex systems metastability occurs, and it requires a time-dependent analysis of the problem at hand, since since considering only the steady-state is no longer sufficient. In many cases this makes analytical treatment very difficult since time-dependent problems are rarely soluble. However, in this case we do have a solution, and in the sections that follow we proceed analytically. In the next section we derive the relaxation time scales to the steady-state for systems exhibiting lock-ins by determining the switching times between the two lock-in states. In the following section we will often refer to decisions as technologies, in line with the technology application drawn above.

Escape from non-optimal decisions
3.2.1 Relating the relaxation time scale to the first passage time. We now look at the lock-in mechanism in more detail using our analytical solution. The general case of metastability in bistable systems is covered by van Kampen [31], who does so in three steps. These steps are vital to the explanation of the lock-in mechanism, and in determining the relaxation time scale l À 1 2 to the steady-state as t ! 1.
Step 1 describes the evolution of the initial distribution. For t ≳ 0, the distribution broadens and probability modes have not yet developed at the equilibria. This is seen in Fig 2(a) for the distribution in red at t = 0.1 (the initial condition at t = 0 being that of equal numbers of agents deciding left and right).
Step 2 is that two distinct modes have developed in P(m(n), t), which decompose into left and right parts respectively. How the distribution has decomposed is dependent on the initial condition and can be seen in Fig 2(a) for the distribution in green. Finally for t � 1, step 3 is that each mode has developed its local equilibrium shape and that the transfer of probability between the modes becomes very small. This is shown in Fig 2(a) for the blue distribution and for the black dots, showing that the distributions between t = 10 and t = 10 10 are indistinguishable. At times t � l À 1 3 , where λ 3 is the third smallest eigenvalue of the master operator in Eq (15), the process can then be effectively modelled as a two-state process between the equilibria, with very small rates of probability transfer between the modes, as illustrated in Fig 2(d). In the following our main interest is in finding these rates of transfer between the behavioural modes and the overall rate of relaxation towards the steady-state.
Before proceeding, we illustrate how to calculate the equilibria values of m, including the unstable one, via the steady-state distribution alone (see Fig 2(b)). For β > β c , the stable equilibria m − and m + are the maxima of the steady-state distribution (−1 and 1 in Fig 2(b)) with the larger of the two modes being the optimal equilibrium value. The intermediate unstable equilibrium point m u is the minima of the distribution (found at N/2 for Fig 2(b)), and a system found either side of m u will likely drift towards the corresponding stable equilibrium point. We will denote the critical number of right deciding agents as n u = N(m u + 1)/2. Note  (16) from an initial condition at n = N/2. Beyond the initial condition we see the emergence of two modes of behaviour, the agents either increasingly choose the left or the right technology. The evolution of the probability distribution becomes very slow for t > 10 and the distribution at t = 10 10 is indistinguishable from that at t = 10. (b) A plot of the analytic steady-state distribution from Eq (22) shows that in the true steady-state limit almost all the agents will be locked into the right technology. Importantly, this is completely distinct from the time-dependent solution even at large times. Note the presence of a small mode at the unfavoured left-hand technology seen in the inset.

PLOS ONE
that for β � β c (see Sec. 4 in S1 File), or for |F| > J(α + 1) in the β � 1 limit (see Sec. 6 in S1 File), there is only one equilibrium value of m.
Defining p L ðtÞ ¼ P m<m u Pðm; tÞ and p R ðtÞ ¼ P m>m u Pðm; tÞ as the probabilities to be in the left and right modes respectively, following the approximative two-state process in Fig 2(d), we can then assert that, where τ lr and τ rl are the mean times to escape the from equilibria at the left and right modes respectively (i.e., starting at the equilibria, how long it takes to reach m u ), and ϕ R is the probability that if the system starts at m(0) = m u that the agents initially coalesce on the right technology. We derive ϕ R in Sec. 5 in S1 File (using similar methods as the calculation of the mean first passage time in the next section). This probability enters the rates of transition to the other mode since τ lr and τ rl only give the times to escape from each mode, which must be scaled by the probability to enter the other mode, since upon reaching m u the agents could also revert back to the decision they inhabited before. Eventually, π L and π R reach their stationary values which are determined by @ t p s L ¼ 0, explicitly, This equation helps us interpret the results of Fig 2(b): for the chosen parameter set τ lr � τ rl : once the agents have made the optimal decision they will not change it on any reasonable time scale. Note that even at the steady-state there is a mode at m = −1 (see inset Fig 2(b)), and that the shape of the left-hand distribution will be near identical to the shape of the distribution at t = 10 in Fig 2(a), but is of much smaller magnitude. The three steps described by van Kampen can also be seen on a SSA trajectory level in Fig 2(c): after an initial period of indecisiveness the individual population trajectories eventually coalesce on either the optimal technology at m = 1 (with slightly greater probability) or else the sub-optimal technology at m = −1.
Before calculating the switching times between the equilibria, one can ask how these switching times are related to the rate of relaxation to the steady-state. Namely, how do τ lr and τ rl relate to the smallest magnitude non-zero eigenvalue λ 2 ? We know that for t � l À 1 3 the only relevant time scale is the relaxation time scale from the metastable state to the steady-state, hence we can write, Pðm; tÞ � P s ðmÞ þ expðÀ l 2 tÞF 2 ðmÞ; ð25Þ where F 2 (m) is the eigenvector of A (see Eq (15)) corresponding to the eigenvalue λ 2 . Inserting Eq (25) into Eq (23) we find that, where we have used the fact that ∑ m F 2 (m) = 0, which comes from summing Eq (25) over all m and enforcing normalisation conditions on P(m, t) and P s (m). As Bouchaud showed in [3], one can approximate the master equation in (14) as a Fokker-Planck equation [31,43], hence finding that the mean first passage times are proportional to t lr / expðNð1 À F=Jð1 þ aÞÞÞ; ð27Þ whose full derivation we show in Sec. 6 in S1 File. The utility of this result is its ease of interpretation: where there are two stable equilibria the mean times for switching between them are exponential in the number of agents. (27) and (28) are very good approximations in the β � 1 and N � 1 limits where m ± = ±1. However, there are some cases where the distributions π L (t) and π R (t) have a greater variance and are less peaked (as illustrated in Fig 3(a)) in which case contributions to the time needed to pass the potential barrier at m u must be taken from multiple values of m. We first define the conditional normalised probability distributions of having m inside each of the the metastable phases as r l ðm < m u Þ ¼ P s ðmÞ= P m<m u P s ðmÞ and r r ðm > m u Þ ¼ P s ðmÞ= P m>m u P s ðmÞ for the left and right equilibria respectively. We determine this from the intuition that the metastable probability modes are the same as those found in the steady-state distribution up to a multiplicative pre-factor, an intuition that we confirm in Fig 3(b). One can then define the mean times to switch between the two equilibrium modes at m < m u and m > m u as weighted sums over the conditional distributions,

A better approximation of the relaxation time scale. Eqs
where τ n is the mean first passage time to reach n u given one starts with n agents deciding on the right-hand technology. In our case it is possible to find the values of τ n exactly, and we do so in line with methods shown in [31,59,60]. Consider again the microscopic transitions previously seen in Fig 1(b). For this birth-death process one can write a backward equation, which is formally the adjoint equation to the master equation, commonly used for first passage processes. For our purposes, it is more convenient to use the discrete time backward equation given by [31,59], where a i = (N − (i − 1))r(i − 1), b i = (i + 1)l(i + 1), Q j,i (t) is the probability of being found with j agents deciding on the right technology a period of time t after being found with i right technology deciding agents, and Δt is the time step (which is taken to zero in the continuous time limit). Note that the absolute time t n at time step n is defined by t n = nΔt. The mean first passage times for hitting m u from m < m u and m > m u must be considered separately. In the following we show the calculation for m < m u , although the procedure is analogous for m > m u . Consider only the states of the system m(n)<m u , where we define m u = m(n u ) with n u being the number of agents deciding for the right technology at the unstable equilibrium point. We define a new Markov process by setting b n u À 1 ¼ 0, which defines m u as an absorbing boundary [18,31]. A schematic of these new dynamics is show in Fig 3(c). From the backward equation for this new process we then define c i ðtÞ ¼ Q n u ;t ðtÞ as the cumulative probability of reaching n u a time t after having i right deciding agents. Hence, the probability of hitting n u at time t is given by ψ i (t) − ψ i (t − Δt), and the mean first passage time τ n is given by a weighted sum over these probabilities, with ψ n (−Δt) = 0. From the backward equation we find the recursive relationship for ψ i (t), which has the intuitive interpretation that the probability of reaching m u at time t + Δt is the probability of hopping to i+ 1 and reaching m u from there, plus the probability of hopping to i−1 and reaching m u from there, plus the probability of not hopping and reaching m u from i [59]. If we subtract ψ i (t) from both sides of Eq (33), multiply by t, and sum over all t we then get a recursion relation for the mean first passage time, One can perform similar calculations for the higher order moments of the first passage time distribution, although we do not consider them here (see [59] for more details). We have two boundary conditions on this recursion relation. The first is t n u ¼ 0, which reflects the fact that if one starts at n u the time taken to reach it is obviously 0. The second is less obvious in that it comes from physical conditions at the left boundary and is τ 1 − τ 0 = −1/a 1 . This reflects the fact that in the state n = 0 there is only one possible move to n = 1, which occurs with average time 1/a 1 . Now, we introduce the difference variable η i = τ i − τ i−1 which transforms Eq (34) into, subject to η 1 = −1/a 1 . This can be solved recursively and the result is, where again we note that the empty product is equal to 1. Finally, to find the τ n from η n we see that, where we have used the fact that t n u ¼ 0. Using the same approach (or via symmetry considerations) one can derive the case of n > n u , whose result is, The results derived in Eqs (37) and (38) are exact. If one wishes for a simpler result, it is possible to simplify the expressions in the large N limit by approximating the product term in Eqs (37) and (38) with ∏ k γ k (n) = exp(∑ k log(γ k (n))) � exp(N R log(γ k (n))dn), where we have explicitly included the dependence on n for clarity (see [60] for more details). Doing so gives similar exponential dependence on N as we explored in Eqs (27) and (28), and hence we do not show the result of the large N approximations on Eqs (37) and (38) here.
We show a plot of τ m(n) � τ n in Fig 3(d) for F > 0, where we clearly see that it is more difficult to escape from the optimal choice technology than the sub-optimal since on average one has to wait much longer to escape the more optimal collective decision. The values of τ n determined via Eqs (37) and (38) can then be used in combination with Eqs (29) and (30) in order to find the relaxation time scale λ 2 in Eq (26). Note that although this determination of λ 2 is still formally an approximation it is very accurate, and we show in Fig 3(e) that for the timedependent solution at t ¼ l À 1 2 it results in a distribution very close to the steady-state distribution. Computationally, using eigvals, we find l À 1 2 � 1288:8 for Fig 3, whereas our approximation gives l À 1 2 � 1279:8, hence only slightly underestimates the relaxation time to reach the steady-state. Our approximation is an underestimation since the relaxation time since it does not take into account the initial relaxation into the metastable state of Oðl À 1 3 Þ. Since for Fig 3 we have N = 50 agents, one can also conclude that our method works well even when not in the large N limit. Moreover, our approximation is even better for systems such as that seen in Fig 2, where the modes at n − and n + are very strongly peaked, since this type of system more closely corresponds to the approximations we have made in our calculation of λ 2 as the time taken to reach the metastable regime is even shorter compared to the time taken to reach the steady-state. The disparity seen between our approximation for λ 2 for the parameters in Fig 3 constitutes a worst case scenario of the calculation of λ 2 . Eq (25) then constitutes an approximate solution for the bimodal regime of the mean-field model that does not require calculation of the fast decaying eigenvalues computationally, valid where t � l À 1 3 . We finally observe that our calculation of λ 2 gives us an approximate time-dependent probability distribution, P(m, t)�P s (m) + exp(−λ 2 t)F 2 (m), which does not require the calculation of any eigenvalues and is computationally valid in the metastable regime where all other eigenvalues of order Oðl À 1 3 Þ or less have decayed away. Practically, in order to implement this result, upon realising that we have λ 1 = 0 and λ 2 we then substitute these into Eq (16) and where t � l À 1 3 we can safely ignore the contributions from i � 3 in the sum.

Model calibration
Brock and Durlauf [1] used their equilibrium solution for econometric analysis and model calibration. However, as we have emphasised in the previous section, an equilibrium solution is often not applicable on any realistic time scale, even for small numbers of agents, and hence time-dependent dynamics of the binary choice model must be considered. Generally speaking, finding an analytic time-dependent likelihood function for an agent-based model is a rarity [25], however our time-dependent solution for the probability distribution allows us to construct one for the mean-field binary decision model, allowing us to conduct practical model calibration in a short time. This makes our solution useful to researchers who do not have access to vast amounts of computing power. We then test our calibration procedure on simulated data with known parameters in order to assess the conditions on the data necessary to provide reliable calibration.

Construction of likelihood function
We begin by defining the likelihood in the standard way. Say we have a set of L data points M ¼ fmðt i Þg for i 2 {1, 2, . . ., L}, measured at times {t 1 , t 2 , . . ., t L } describing the evolution of the order parameter m(n(t)) over some time period. In general one cannot assume this process is at equilibrium and the evolution from the initial condition at m(t 1 ) does not follow a steadystate trajectory. This requires us to use Eq (16) in order to calculate the probability of observing the trajectory M given some assumed parameter set θ = {F, J, γ}, known as the likelihood, and is given by, where P θ (m(n(t i ))) is Eq (16) evaluated for parameter set θ at time t i . The likelihood function is expressed as a product of probabilities over the whole time series since it is the probability of observing m(t 1 ) at t 1 and m(t 2 ) at t 2 and so on (hence from the laws of probability is multiplicative). Note that since β pre-multiplies F and J and (1 + α) pre-multiplies βJ in the utility function, β or α cannot be inferred separately but only βF and β(1 + α)J can be inferred [1]. Hence for calibration we set β = 1 and α = 0 and infer only F, J and γ. The aim of the calibration procedure is then to find the parameter set θ ? that maximises LðMjyÞ with respect to θ. Generally, since the likelihoods are generally very small quantities, it is more computationally convenient to instead minimise the negative log likelihood À lnðLðMjyÞÞ which generally takes values �1. In most cases (including ours), the parameter set corresponding to the minimum value of À lnðLðMjyÞÞ cannot be analytically determined and hence one must proceed algorithmically by testing many different values of θ via an optimisation algorithm. In this paper we utilise the adaptive differential evolution optimiser from the Julia package BlackBoxOptim [61,62], which is determined to be the best likelihood optimisation algorithm compared to several other state-of-the-art methods [62]. Using this algorithm, the optimal parameter set θ ? is then determined through, where Θ is the set of all possible parameters in the optimisation range selected. In the calibration that follows the true parameter set for the simulated data is θ true = {F = 0.025, J = 1.5, γ = 1.0} (same as in Fig 3) and the parameter ranges of optimisation for the set Θ were chosen such that F 2 [−2, 2], J 2 [exp(−2), exp (2)] and γ 2 [exp(−1), exp(1)], where parameters J and γ are optimised in log-space. For the purpose of analysing the errors in the calibration procedure below we define, where E tot is the total error on the inference procedure and f is the error in the fraction between F and J. In most econometric analyses it will be the error in f which is most important since it tells us whether one has properly inferred the magnitude of the influence of exogenous versus endogenous effects on the agents. These error functions will then allow us to determine the conditions necessary of data such that one can conduct a valid model calibration. We note that θ ? does not generally correspond the global minimum of À lnðLðMjyÞÞ, but instead a good local minimum. This is due to the complexity of the function À lnðLðMjyÞÞ, a problem for which there is no known solution to produce the global optimum in a reasonable computational time [63]. The number of optimisation steps we choose for the adaptive differential algorithm is 500, with a population size of 1000 such that each calibration procedure took less than 1500s on a laptop with Intel Core i7-8650U CPU 1.90GHz x 8 running Ubuntu 18.04.5 LTS.

Calibration from a single trajectory
In many real-world situations data is only collected once for a particular event. For example, national elections do not have multiple realisations, hence if one was to collect the data of the percentage of votes for Democrats or Republicans at all elections from 1868 (from when only Democrats or Republicans have won the presidential vote) then one ends up with a single trajectory of data. Therefore, it is important to assess the conditions under which calibration procedures are accurate with respect to data for which the parameters are known. For this purpose we utilise the SSA to provide stochastic trajectories for a particular realisation of our mean-field system for the parameter set explored in Fig 3, which notably expresses bimodality. First, we look at the case of calibration where only one of the two equilibria is explored, from the trajectory seen in Fig 4i(a). The question here is: is there enough information in a single trajectory (from a bimodal system) exploring a single equilibrium mode for the correct parameters to be inferred? The answer is no, which can be seen from a plot of the errors in the inference against the calibration time (for a fixed number of 100 data points) used for the inference in Fig 4i(b) and 4i(c). Even though increasing the calibration time improves the calibration, the overall inference of the parameters is still poor for large calibration times with the errors being � 1. The reason for this is simple: for a system that expresses bimodality the probability is shared across both modes, hence for a trajectory that only realises one of its two equilibria in a given realisation, the optimiser will instead choose a set of parameters θ ? that corresponds to the single mode explored in the data.
Having ruled out the ability to conduct good optimisation when only a single equilibria is explored, we now look at the case where a system realises both of its modes of behaviour on a given trajectory, shown in Fig 4ii. The trajectory for this realisation is shown in Fig 4ii(a), where the mode of behaviour changes approximately half way through the trajectory. As seen from Fig 4ii(b) and 4ii(c), the calibration is significantly improved when both equilibria are explored in the calibration procedure (using a fixed number of 100 data points as before). The reason for this is relatively clear: parameter sets for which the probability distribution exhibits only a single mode near m = −1 or 1 have a lower likelihood, when the system explores the opposite mode from the one they describe. Hence, one can conclude that if the system is bimodal, and the data expresses both of these behaviours, the calibration produced will be relatively good.
Finally, one may want to know the optimal number of data/time points to use, rather than the optimal time of calibration. Fig 4iii(a) shows the cases of calibration using differing numbers of data points over the entire trajectory considered in Fig 4ii(a), from 1001 to 10 data points. Fig 4iii(b) and 4ii(c) show the results of inferring the parameters with respect to the errors E tot and f: one observes that having more data points is better, but recognises diminishing returns with the additional data points beyond around 10 2 . Therefore, having more data points generally leads to a better parameter estimate but is not as important as having data that represents all behavioural modes of a given parameter set.

Calibration from multiple trajectories
On some occasions there exist socio-economic binary decisions with more than one realisation. For example, although American presidential elections only happen once every four years, prior to the elections pollsters survey the public to assess what the current opinion is [64]. These multiple surveys across the different elections constitute different independent realisations of the same socio-economic phenomena, if one excludes pollster bias or the possibility that over different elections the underlying model parameters are different. Hence, having multiple realisations of the same phenomena does not mean the existence of multiple realities, but simply temporally or spatially separated independent events which one assumes have similar exogenous and endogenous influences on the agents. Having different realisations can be very beneficial for the calibration procedure, since as we saw for the case of a single trajectory, if only one behavioural mode is explored the calibration is weighted towards parameter sets that

PLOS ONE
favour that reality. Note that the example we explore the initial condition is fixed and the independent trajectories have exactly the same underlying parameter sets (in reality there would be some variations in F, J and γ between the independent events). However, it shows the considerable effect that multiple realisations have on model calibration.
In Fig 5(a) we show 1000 realisations of the process for the parameters in Figs 3 and 4 (in yellow), highlighting five of the trajectories (in black). Fig 5(b) and 5(c) show that for an increased number of realisations the error in the calibration is much reduced, for both E tot and f. This reduction in error is even more pronounced than the reduction due to increased calibration times or number of data points seen in Fig 4. The reason for this is that if one has multiple trajectories of the same process then it is much more likely that both equilibria will be explored, and also that the number of trajectories that explore each equilibria is proportional to the probability that a single trajectory will be found in that state.
In summary, we find that there are several conditions on data that can result in more accurate model calibration, which are: 1. If the system can express bimodal behaviour, and one only has a single data trajectory, then it is imperative that the trajectory used explores both of these behavioural modes. If this is not the case the calibration procedure will weight itself towards parameter sets that only realise the behaviour seen in the data.
2. Having a longer calibration time, or an increased number of data points, both improve the accuracy of the calibration, but with diminishing returns.
3. Having access to more than one realisation of the situation one wishes to model is highly beneficial and results in very accurate calibration. This is seen in Fig 5(b) and 5(c) where the number of trajectories exceed approximately 10 2 . Note also the sizeable drop in error (in log-scale) from a data-set with one realisation to a data-set with 11 realisations in both of these figures. We stress that from a single realisation of real-world data one does not know a priori whether the system expresses bimodality, hence having multiple realisations can be very informative, and provide for much better calibration.

Other methods for model calibration
Although it is very convenient for us to utilise the likelihood function based on our timedependent solution, other likelihood-free calibration methods are available. One of particular mention is approximate Bayesian computation (ABC), which generates SSA (Monte Carlo) trajectories for a given parameter set θ which can then be compared to the original data using a distance function [65,66]. The choice of distance function is important, and typical examples include the Euclidean distance between the trajectories, the Hellinger distance, or even the Wasserstein distance [67]. Often it is best to experiment with each of these, or combinations of them, to see empirically which gives the best result from simulated data, before applying the calibration procedure to a real-world data-set.

Neoclassical economics versus the binary decision model
In this section we compare the results of the paper to the tenets of neoclassical economics (NCE), and in particular explore: (i) information and expectations, (ii) the dichotomy between the social planner and altruistic agents, and (iii) the implications of our findings on model calibration to economics. Note that many of the points we make here are not first noted by us, and one may additionally want to read further publications such as [3,6,37,39,[68][69][70] which have inspired this author.

Information and expectations
NCE assumes that agents have perfect rationality in the sense that all agents can instantaneously perform simultaneous complex calculations and arrive immediately at an equilibrium based on a total knowledge of present and future [39]. However, in the binary decision model no two agents make a decision at the same time, and the maximum foresight an agent can have is to choose the decision that maximises the agent's utility at that point in time. There is no capacity for the model agents to look into the future, or else to gain an insight into which technology is the optimal one, and hence lock-ins still occur within the binary decision model even where agents are 'perfect' (in the limit β ! 1). This would not occur in a neoclassical model as the agents would all decide to choose the optimal technology in order to maximise their individual utilities.

Social planner versus altruistic agents
In NCE the ideas of having selfish agents and a social-planner/auctioneer go hand-in-hand. However, in the binary decision model we have explored, the 'social planner' situation more closely corresponds to a system of altruistic agents of altruistic strength α = 1 rather than the system of selfish agents, since the altruistic agents at least consider how each agent can change to best affect the global utility. This dichotomy is curious, since from the perspective of our agent based model the social planner and selfish agents cannot co-exist. As Bouchaud noted in [3], the papers of Grauwin et al. [42,71] are remarkable in that they exhibit situations in which selfish agents fail spectacularly in a Schelling-like global coordination problem, breaking Adam Smith's 'invisible hand', whereas having agents that are even somewhat altruistic remedies the situation for the benefit of all agents. One can further ask a similar question for the mean-field binary decision model: does having altruistic agents increase the ability of the population to increase global utility? The answer to this question is mixed. In one sense the answer is yes-increasing α increases the endogenous influence and hence the time taken to coalesce on a single decision in the population is much reduced. Additionally, in each behavioural mode the fluctuations are reduced if agents are altruistic, shown in the right-hand plot in Fig 6, decreasing global utility compared to selfish agents. However, in another sense the answer is no-altruistic agents do not know a priori which technology, left or right, is the optimal one and hence can easily coalesce on the wrong technology. In the situation where this occurs the waiting time for them to leave the non-optimal decision is much greater than the selfish agents (shown on the left of Fig 6), and in fact the waiting time to equilibrium is roughly exponential in α. Hence as Brock and Durlauf [1] state, altruistic agents are more susceptible to conformity effects; however, their ability to reach the optimal equilibrium mode often does not occur on a relevant time scale.

Implications of calibration procedure
Previously in Section 4.3, we defined the conditions for good model calibration, which in essence were: many data points, long calibration times, and multiple trajectories. However, the opinion of some of the literature [25] seems to be that one of the main issues with model calibration lies in not having explicit likelihood functions. In our case, even with an analytically defined likelihood function, and only three parameters to infer, model calibration is still a difficult task unless one utilises a well-informed data-set. Additionally, as we have already commented, it is possible to conduct calibration without an explicit likelihood function using likelihood-free methods for parameter inference [66]. It is hence vital for model calibration that either calibration methods become much smarter, for example using methods beyond simply using a likelihood function as done in [72] (here in the context of gene expression), or else the data used for calibration becomes more informative. Although it was found that multiple data realisations are the biggest improvement that can be made to the inference procedure, in practise, for many questions of interest, it is not possible to access such data. Hence, the The left-hand plot shows that as agent altruism increases the metastable waiting times (black line) become exponentially large in α. The blue dots show an exponential approximation to l À 1 2 , clearly exhibiting its close-to exponential behaviour. The right-hand plot shows the change in the equilibrium (t ! 1) distributions. For increasing agent altruism populations become more polarised in each behavioural mode (with much reduced fluctuations) even in the metastable regime.
https://doi.org/10.1371/journal.pone.0267083.g006 need for quantities of good quality economic data is a pressing one since economists are already using such calibrated models for economic prediction. In lieu of such well informed data-sets with which to conduct econometric analysis, economists are limited to optimising their models for one of many indistinguishable parameter sets from the calibration of their models.

Conclusion
In this paper we have solved the binary decision model of Brock and Durlauf [1] in time by mapping it to a stochastic birth-death process which can be solved via the method of [44]. We explored the waiting times to the steady-state in the metastable regime using a very accurate method based on first passage time theory, and used the solution to construct the likelihood function for use in model calibration. The method we employed from [44] can also be used to solve Kirman's ant recruitment model in time [7,10,11]. We additionally note that our solution in Eq (16) is not only a solution to the binary decision model, but is also a time-dependent solution to the mean-field Ising model used as an approximate description of magnetic behaviour in physics. Previous time-dependent solutions to this problem have only come in the form of hydrodynamic approximations [73]. One of the main utilities of our solution is the ability to infer {F, J, γ} from real economic binary decision data, and that these three parameters completely specify the dynamics of the agents (up to the choice in decision rule and multiplicative factors β and α). We constructed a likelihood function and showed that model calibration, even in our simple mean-field model, is a non-trivial task unless one has access to an informative data-set made up of multiple realisations of the same socio-economic phenomena.
Our solution provides a platform for several enhancements and investigations of the model, including: 1. An exploration of non-logit decision rules (such as the Arrhenius form), and additionally decision rules that break detailed balance. This can easily be done using our solution in Eq (16), since it is general in the form of the transition rates W n (S i ! −S i ). The exploration of non-detailed balance decision rules is of particular importance, as emphasised by Bouchaud [3], since economies are not closed systems but are subject to energy influx/out flow. Recent publications indicate that breaking detailed balance results in system properties that are not consistent with properties of systems pertaining to detailed balance [74].
2. In this paper we have not explored our solution with respect to time-dependent F(t) (whose solution is discussed in Sec. 2 in S1 File), in particular with the rise of metastable states. An understanding of how the time-dependence on F(t) can affect metastability and lock-in effects would have relevance to policy makers, for example in affecting the transition to more renewable energy sources by reducing the time spent fixated in the non-optimal side of a technology lock-in. A question of interest would be: given a system of agents is in the lock-in state in the non-optimal technology, what is the optimal time-dependent function F (t) that releases the agents from this lock-in (in minimal time) whilst keeping R 1 0 FðtÞdt to a minimum (i.e., reducing the effort necessary to break the lock-in)?
3. Developing methods to analytically solve multiple choice models, such as those explored in [2,34]. These models are more generalised than the binary decision model we have discussed in this paper, and in fact contain the binary decision model as a special case, so solutions for multiple choice models would have more potential applications. 4. A more detailed investigation of model calibration using more realistic simulated data-sets and real world data-sets. Can one do better than the classic likelihood based inference using more advanced methods? Additionally one could explore the performance of likelihoodfree methods of calibration on models for which we do not yet have analytic solutions.
In our opinion, it is the final item in the above list that is of most importance for future investigation, since as we have found even for three parameter inference of {F, J, γ} in our model, calibration is a difficult task without suitable data-sets. The collection of such data-sets (beyond that of the stock market) that can be used for calibration is therefore of great importance, especially given the need for policy makers to understand immediate challenges relating to the climate crisis and the transition to a more sustainable world [75]. As stated by Beinhocker in The Origin of Wealth [39], economists will get no sympathy from biologists and physicists who must go to great lengths to collect data to test their theories. Socio-economic modelling, and the models similar to that studied in this paper, can form an essential part of the work necessary to understand what policy makers need to do to change human collective behaviour in uncertain times.